#############################################################
## Description: combine sdm and range rasters
##
## Date: 2018-11-17
## Author: Fabian Fopp (fabian.fopp@usys.ethz.ch)
#############################################################


combine.range.sdm.fun <- function(species_name, rast_sdm, rast_range, nmodel=6, final_resolution, dest_plot, dest_raster, write_plot=F, return_raster=T, write_raster=T, overwrite=F){
  
  library(raster)
  
  if(class(rast_range)!='RasterLayer' | class(rast_sdm)!='RasterLayer'){
    warning('Range or sdm raster does not exist. No output file created')
    return(NULL)
  } else {
    
    if(final_resolution<res(rast_range)[1] | final_resolution<res(rast_sdm)[1]){
      warning('Resolution of range raster (', round(res(rast_range)[1],3), ') or sdm raster (', round(res(rast_sdm)[1],3), ') is coarser than desidered final resolution (', final_resolution, ').
                ', round(order(final_resolution<res(rast_range)[1], final_resolution<res(rast_sdm)[1])[2],3), ' will be used as final resolution')
      final_resolution <-  res(rast_sdm)[1] #is rast_range really completely irrelevant?
    }
    
    #rast_range <- projectRaster(rast_range, rast_sdm) #gives weird results (range raster cells will not be big squares anymore but can have circular shapes)
    #origin(rast_sdm) <- origin(rast_range)
    #rast_range <- disaggregate(rast_range, fact=res(rast_range)/res(rast_sdm)) #also not perfect, but close
    #rast_range <- projectRaster(rast_range, rast_sdm) #now I have values between 0 and 1!
    rast_range <- resample(rast_range, rast_sdm)
    rast_range <- rast_range >= 0.5
    rast_sdm <- rast_sdm >= nmodel
    rast_combined <- rast_range + rast_sdm
    rast_combined <- rast_combined >= 2

    if(res(rast_combined)[1]<final_resolution){
      rast_combined <- aggregate(rast_combined, fact=final_resolution/res(rast_combined)[1])
      rast_combined <- rast_combined>=0.5
    }
    
    
    #plot combined map
    if(write_plot==T){
      if(!dir.exists(dest_plot)){
        dir.create(dest_plot, recursive=T)
      }
      if(file.exists(file.path(dest_plot, paste0(species_name, ".jpg"))) & overwrite==F){
        stop('Plot already exists. Use overwrite=T to overwrite it.')
      }
      jpeg(file.path(dest_plot, paste0(species_name, ".jpg")), height=2000, width=3000)
      plot(rast_combined, col=c(rgb(0,0,0,0), rgb(0,1,0,1)), legend=F, main='range and sdm combined', cex.main=3, axes=F, box=F) 
      map('world', add=T)
      dev.off()
    }
    
    #save species range raster
    if(write_raster==T){
      if(!dir.exists(dest_raster)){
        dir.create(dest_raster, recursive=T)
      }
      if(file.exists(file.path(dest_raster, paste0(species_name, ".asc"))) & overwrite==F){
        stop('Raster already exists. Use overwrite=T to overwrite it.')
      }
      writeRaster(rast_combined,filename=file.path(dest_raster, paste0(species_name, ".tif")), overwrite=TRUE, datatype='INT1U')
    }
    
    if(return_raster==T){
      return(rast_combined)
    }
    
    rm(rast_combined, rast_range, rast_sdm)
    cat("########### End of computation for species: ",species_name," #######################", "\n") 
  
    gc()
  
  
  }
}